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ABSTRACT 


This  research  study  is  concerned  with  spatial  (relative) 
phase  reconstruction  using  an  imaging  sensor.  Specifically, 
work  was  completed  which  investigates  the  retrieval  of  the 
spatial  phase  of  the  aperture  wavefront  from  sampled  inten¬ 
sity  data  in  the  focal  plane  of  an  imaging  sensor,  where  phase 
is  not  explicitly  present. 

Techniques  investigated  are  not  interferometric  methods, 
hut  are  signal  processing  techniques  that  use  the  Fourier 
transform  properties  of  a  lens,  analytics,  and  numerics  to 
retrieve  the  spatial  phase  of  the  aperture  field  from  inten¬ 
sity  measurements  in  the  focal  plane  of  an  imaging  sensor. 

The  techniques  investigated  include  an  analytic  approach 
which  is  used  to  develop  several  numerical  approaches,  and 
the  Gerchherg-Saxton  Algorithm  which  is  adapted  to  solve 
this  specific  phase  retrieval  problem.  Finally,  the  results 
from  simulations  of  these  methods  are  compared  yielding  the 
Gerchherg-Saxton  Algorithm  as  the  most  promising  method  to 
approximate  the  relative  phase  in  the  aperture  plane. 
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PHASE  RETRIEVAL  USING  AN 
IMAGING  SENSOR 

I  Introduction 

The  performance  of  state-of-the-art  imaging  sensors 
is  limited  by  anomalies  such  as  the  turbulent  atmosphere 
and  aberations  introduced  by  the  optics  of  the  sensor.  As 
optical  radiation  propagates  outward  from  its  source  the 
turbulence  (random  fluctuations  in  the  refractive  index)  in 
the  atmosphere  distort  the  shape  of  the  optical  wavefront 
as  well  as  cause  intensity  fluctuations  across  the  optical 
wavefront.  The  random  distortions  in  the  shape  of  the 
optical  wavefront  are  the  primary  contributors  to  the 
degradation  of  image  quality  in  imaging  systems  (Refs  1  and 
2). 

In  order  to  maximize  image  quality  and  resolving 
power  of  an  imaging  sensor,  real-time  active  optical  systems 
can  be  used  to  compensate  for  the  distortions  in  the  shape 
of  the  optical  wavefront,  adjusting  the  optical  path  length 
of  the  rays  of  each  wavefront  to  be  equal.  Active  optical 
systems  consist  of  three  main  components?  these  are«  a 
wavefront  sensor  (estimates  the  distortion  in  the  optical 
wavefront),  active  optics  (such  as  deformable  mirrors)  in 
the  imaging  sensor,  and  a  feedback  loop  (Ref  2).  A  Wave- 
front  sensor  is  an  important  component  in  closed-loop 
active  optical  systems.  This  research  project  is  concerned 
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with  using  signal  processing  techniques  to  approximate  the 
optical  wavefront  distortions  incident  on  an  imaging  sen¬ 
sor. 

Background 

State-of-the-art  imaging  systems  use  interferometers 
to  sense  wavefront  distortions,  since  detectors  measure 
intensity  of  the  electric  field  and  can  not  directly 
measure  the  phase  of  a  field.  The  basic  configuration  of 
an  imaging  system  using  an  interferometer  to  sense  wave- 
front  distortions  is  shown  in  Figure  1.  The  1:  itations  of 
these  active  optical  systems  are  that  they  require  two  sen¬ 
sors,  one  for  wavefront  estimation  and  a  second  sensor  for 
imaging,  adding  complexity  to  the  system  increasing  the 
size,  weight,  and  cost  (acquisition  and  maintenance).  The 
use  of  a  beamsplitter  (inherent  in  this  class  of  imaging 
systems)  makes  the  system  more  susceptible  to  noise  (Refs  1 
and  3)  ■ 

Currently  there  is  interest  in  developing  a  signal 
processing  technique  that  can  be  used  to  retrieve  the  phase 
of  the  aperture  field  from  intensity  measurements  in  the 
focal  plane  of  an  imaging  sensor,  where  phase  is  not  explic¬ 
itly  present  (Fig  2).  There  are  two  distinct  advantages 
for  using  signal  processing  techniques  instead  of  inter¬ 
ferometric  techniques  for  wavefront  sensing  in  imaging  sen¬ 
sors.  The  first  advantage  is  that  only  one  sensor  is  re¬ 
quired  for  imaging,  since  intensity  data  in  the  focal  plane 
can  be  measured  with  state-of-the-art  detectors.  The  second 
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Figure  1,  Basic  configuration  of  an  imaging  system  using  an 
interferometer  for  wavefront  sensing 
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Figure  2.  Basic  configuration  of  an  imaging  system  using 
signal  processing  of  intensity  data  to  retrieve  the  phase 
of  the  wavefront 


advantage  is  that  the  signal  processing  technique  does  not 
require  the  use  of  a  "beamsplitter,  and  does  not  reduce  image 
intensity  which  is  especially  important  in  systems  where  the 
intensity  of  the  source  can  not  be  increased  or  can  only  be 
increased  at  the  expense  of  the  flexibility  of  a  system 
(size,  weight,  and  cost). 

The  problem  of  estimating  the  phase  of  an  aperture 
field  from  intensity  measurements  in  the  focal  plane  is  not 
a  new  one  and  is  commonly  referred  to  in  the  literature  as 
the  phase  retrieval  problem.  This  phase  retrieval  problem 
has  many  applications  in  other  areas,  such  as«  x-ray  crys- 
allography,  particle  scattering,  optical  communications, 
radar,  and  lens  testing  (Refs  2, 4, 5. 6, 7, 8,  and  9). 

Problem  Statement 

The  object  of  this  thesis  is  to  study  techniques 
which  can  be  used  to  retrieve  the  spatial  (relative)  phase 
of  a  complex  (consisting  of  real  and  imaginary  parts)  wave- 
front  incident  on  an  imaging  sensor,  from  intensity  measure¬ 
ments  in  the  focal  plane  of  this  sensor. 

Approach 

This  research  report  is  organized  into  five  chapters i 
the  Introduction,  Formulation  of  the  Problem,  Phase  Retrieval 
Techniques,  Simulation  Results  and  Analysis,  and  the  Conclu¬ 
sion.  The  second  chapter,  Formulation  of  the  Problem, 
develops  the  theory  needed  to  retrieve  the  phase  of  the 
aperture  field  from  intensity  measurements  in  the  focal 
plane  of  an  imaging  system.  The  third  chapter,  Phase 


Retrieval  Techniques,  presents  the  two  approaches  used  in 
this  research  study  to  infer  the  phase  of  the  aperture  field; 
these  are;  The  Phase  Analytic  and  Numeric  Techniques,  and 
the  Gerchberg-Saxton  Algorithm  with  and  without  feedback. 

The  fourth  chapter,  Simulation  Results  and  Analysis,  defines 
the  aperture  wavefronts  used  in  the  simulation,  the  criterion 
used  to  rate  the  performance  of  the  phase  retrieval  techniques 
simulated,  and  finally,  presents  the  results  and  analysis  of 
these  phase  retrieval  techniques  simulated  as  a  function  of 
one  spatial  variable.  The  fifth  chapter,  the  Conclusion, 
is  a  summarization  of  the  results  of  the  research  accom¬ 
plished  during  this  research  study  and  recommends  areas  for 
further  investigation. 


II  Formulation  of  the  Problem 
The  purpose  of  this  chapter  is  to  develop  the  theory 
needed  to  retrieve  the  phase  of  the  aperture  field  from 
intensity  measurements  in  the  focal  plane  of  an  imaging 
system.  This  research  is  concerned  with  the  retrieval  of 
the  phase  of  a  wavefront  (U  (xo ,y„ ) ) ,  incident  on  the 

a  a  a 

aperture  of  an  imaging  sensor,  from  intensity  measurements 
(Im(x^,yf) )  in  the  focal  plane  of  this  sensor. 

It  is  known  from  Fourier  optics  that  the  relationship 
between  the  focal  plane  (U^.(x^,yf))  and  the  aperture  plane 
(U  (x  ,y  ))  electric  field  (Fig  3)  is  given  by 

a  a  a 

*  (jTf) 

<V  ' 

ff  (l) 

-Cp  ***  ' 

The  spatial  Fourier  transform  is  defined  to  be 

c  JfU.Ky^( if"  ¥*)) (2) 

Substituting  equation  (2)  into  equation  (1),  the  focal 
plane  field  is 

IfrU,*)*  (p)  filled)  (3) 
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stem 


where 


\  =  the  wavelength  of  the  radiation 
f  =  the  focal  length  of  the  imaging  sensor 
f  =  spatial  frequency  in  the  focal  plane  (x  direction) 
f  =  spatial  frequency  in  the  focal  plane  (y  direction) 

v 


Equation  (3)  is  valid  for  monochromatic  (single  wavelength) 
radiation  as  well  as  quasi-monochromatic  radiation.  The 
conditions  that  must  be  met  for  radiation  to  be  considered 
quasi-monochromatic  are*  the  bandwidth  (BW)  of  the  radi¬ 
ation  must  be  much  less  than  the  center  frequency  (vQ)  of 
the  radiation 


BV/ 

~v7 


<  <  1 


(4) 


and  the  reciprocal  of  the  bandwidth  must  be  greater  than 
the  length  of  the  longest  optical  path  (d)  involved  divided 
by  the  speed  of  light  (c) 


BW 


(5) 


(Ref  9). 

We  are  concerned  with  an  imaging  sensor  so  intensity 
data  in  the  focal  plane  (lm(xf,yf))  is  a  measured  quantity, 
and  is  the  modulus  of  the  electric  field  (U^x^ty^)) 


=  I  Wf  r  Vf(Xt,h)U*(*4.'f*) 


(6) 


where  the  asterisk  (*)  is  used  to  denote  the  complex  conju¬ 
gate.  The  phase  terms 


(j )  exp  l^sr)  <?xp  (xf+yS)) 


in  equation  (3)  can  be  dropped  for  convenience  since  they  are 
of  no  consequence  in  the  measurement  of  the  intensity  in  the 
focal  plane.  For  simplicity  the  relationship  between  the 
focal  plane  and  the  aperture  plane  electric  fields  can 
be  defined  by  the  spatial  Fourier  transform  relationship 
(except  for  a  constant  term),  specifically 

<*>  OP 

Uf  (Mf)  -  (z?)  ffthfay 


-OO  -CO 


*(»)  li 


to  CO 

=  (a?)  ?*/  [uffa, Yf)j 


For  a  more  thorough  explanation  of  Fourier  optics  and  the 
Fourier  transform  properties  of  a  lens,  refer  to i 
“Introduction  to  Fourier  Optics",  by  J.W.  Goodman  (Ref  9)* 
For  mathematical  simplicity  it  is  assumed  that  the 
aperture  wavefront  is  of  unit  magnitude  and  has  a  phase 
component  (wo (x„ , y  ) ) ,  where 

cl  o.  el 

UA**>1*)  '  P^y-j  (12 

The  pupil  function  (P(xQ  y  ))  is 

d- 1  a. 


p 


j1  <?/se 


(13) 


where 


wx  =  width  of  rectangular  aperture  in  x  direction 
w  =  width  of  rectangular  aperture  in  y  direction 

v 

The  phase  component  (w  (x  ,y  ))  is  the  phase  distortion  in 

a  cl  a 

the  aperture  field  that  is  to  be  retrieved  in  this  research 
study.  Now,  substituting  the  aperture  wavefront  (equation 
(12))  into  equation  (8)  an  expression  for  the  focal  plane 
field  is  found 
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co  eo 

(if) 

~  -  co  -’at 

eXP  ("-tr  (X<X<x+  ^  V°-)) dx^dya.  (14) 

Using  equation  (6)  the  intensity  (Im(xf,yf))  in  the  focal 
plane  is  found 


X*  (**  >  /f ) =  fyjj'  ft  k  ►)  ft- «  )  etfifo  Mi**  /«•)-  Uk«  /*.' )) 

*•>) * /f  fr0'" y* j)i AfiixUyl (15) 


Equation  (15)  illustrates  that  even  though  there  is  no 
explicit  phase  information  in  the  intensity  measurement 
in  the  focal  plane  of  an  imaging  sensor  this  intensity 
distribution  is  non- linearly  related  to  the  phase  of  the 
aperture  field. 

The  development  of  the  theory  needed  to  investigate 
the  retrieval  of  the  phase  of  the  aperture  field  from  inten¬ 
sity  measurements  in  the  focal  plane of  an  imaging  system 
is  now  complete.  The  ideas  and  theory  developed  in  this 
chapter  are  used  in  the  phase  retrieval  techniques 
investigated  in  this  thesis. 
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Ill  Phase  Retrieval  Techniques 


The  purpose  of  this  chapter  is  to  present  the  phase 
retrieval  techniques  investigated  in  this  research  project. 
The  chapter  is  divided  into  two  sections.  The  first  section 
develops  the  Phase  Analytic  and  Numeric  Techniques.  These 
techniques  represent  the  aperture  plane  wavefront  in  a 
finite  Fourier  series  and  the  approximates  the  Fourier 
series  coefficients  as  a  method  to  reconstruct  the  aperture 
field.  The  second  section  develops  the  Gerchberg-Saxton 
Algorithm  and  a  modified  input-output  approach.  These 
Gerchberg-Saxton  Algorithm  techniques  directly  retrieve 
the  aperture  phase,  i.e.  they  do  not  compute  the  Fourier 
series  coefficients. 

Phase  Analytic  and  Numeric  Techniques 

Analytic  Approach.  The  first  approach  is  directed 
at  developing  an  analytical  solution  to  the  phase  retrieval 
problem  and  is  referred  to  in  this  paper  as  the  Phase 
Analytic  Technique.  The  first  step  in  the  Phase  Analytic 
Technique  is  to  represent  the  aperture  wavefront  in  a 
finite  Fourier  series.  Using  a  known  pupil  function 
(equation  (13))  and  the  representation  for  the  focal  plane 
wavefront  (equation  (12))  an  expression  for  the  aperture 
wavefront  is  obtained  (Ref  9) .  Now,  expanding  the  aperture 
field  in  a  Fourier  series 


(16) 


P(^/a)ex^ G  ~  Pt**^*)  Fn/a 

exp  (^  z7T((  )X**fwyV*)  J 

where  the  (F  )  terms  are  complex  Fourier  coefficients. 
After  substituting  the  above  expansion  of  the  aperture 
field  (equation  (l6))  into  equation  (8)  an  expression  for 
the  focal  plane  field  is  obtained 


After  performing  the  above  transform  we  find  the  focal 
wavefront  to  be 


For  notational  simplicity  the  sine  function  is  defined  to 
be 


S.MC  X  ^ 


■&r‘ (*«•#)) 


5/A/ 

_ _ 

TTW*  tx  f  s£f) 


(19) 


Using  this  sine  function  (equation  (19))  equation  (18) 
simplifies  to 
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CD  OD 

(4(*,*) s  ir*Z  Z  £*  s,NC«  &  swc«  y#  (20) 

N--00  /YU-* 

Equation  (20)  forms  an  orthogonal  set  of  basis  functions 
to  represent  the  focal  plane  field.  From  equation  (20)  an 
expression  for  the  focal  plane  intensity  is  found 


Ifiw*)  *  I  jrZ Z ^  s'^h  s"*=«yf  I  (21) 


fts-to  A\»*® 


Now,  we  have  a  representation  for  the  focal  plane 

intensity  (If(x^,yf))  and  we  know  the  exact  focal  plane 

intensity  (Im(xf ,y^ ) ) ,  which  is  a  measured  quantity.  For 

optimization  purposes  a  method  for  determining  the  Fourier 

series  coefficients  (F  )  is  needed.  A  cost  function  (J) 

nm 

is  defined  to  minimize  the  mean  square  error  between  the 
focal  plane  intensity  representation  (lf(x^.,y;f))  and  the 
measured  focal  plane  intensity  (Im(xf ,yf ) ) ,  specifically 

<b  CO 

I )  ~-£f  fa* ft )  I  dxfdyj  ^22^ 

-OJ  “OO 


This  cost  function  (J)  is  known  as  the  "approximation  by 
the  method  of  least  squares  or  an  approximation  in  the 
mean”  (Ref  10).  Expanding  equation  (22)  yields 


J-jj  lljtf.it) 


15 


T 
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After  substituting  the  representation  for  the  focal  plane 
intensity  (equation  (21))  into  this  cost  function 
(equation  (23)),  a  closed  form  solution  to  this  optiza- 
tion  problem  could  not  be  found.  An  excessive  number  of 
cross  terms  are  generated  by  the  representation  of  the 
intensity  in  the  focal  plane  squared  ( (I^(xf ,yf ) )2)  term. 
Specifically 

■If  (*»•&)*  sinc„yf 

£. fa su,t* I*  <24) 

pa-a>  Qa  -CD 

The  fact  that  the  representation  of  the  focal  plane  electric 
field  is  an  orthogonal  series  does  not  help  in  this  case 
since  we  are  dealing  with  powers  of  intensity. 

Numerical  Approaches.  An  analytical  closed  form 
solution  could  not  be  found,  so  two  numerical  approaches 
were  developed.  These  numerical  approaches  use  the  measured 
intensity  in  the  focal  plane  to  generate  the  Fourier  series 
coefficients.  These  estimated  Fourier  series  coefficients 
are  then  used  to  generate  an  approximated  complex  aperture 
wavefront. 

The  first  numerical  approach  is  the  Phase  Numeric  1 
Technique.  This  technique  uses  equation  (20)  to  find  an 
expression  for  the  Fourier  series  coefficients,  specifi¬ 
cally 


00  CO 


J  JUf  (Xf,Y<)  5/f/C*  ft  dx;  Jy,  - 

5//^/,  S/*C*ft  S/A/c24^crn4Jy<(25) 


U)x 

-  00  -So  N.  -  00  M*  -OS 

Solving  equation  (25)  for  the  Fourier  series  coefficients 
yields 

a o  ao 

I  ,  .  .  V  _  .  *  ..  _  K  ..  /  / 

(26) 


ao  ao 

f~xj-  -  (\f^  ^ ^  (J-f  Obi  Yf)  SlNC*  Yf  S/MCj  Kf  c/xfdy^ 


-CO  -CO 


which  is  an  exact  expression  for  the  Fourier  series 
coefficients.  An  approximation  is  made  by  using  the 
modulus  of  the  focal  plane  field  (from  measured  intensity 
data)  instead  of  the  actual  unknown  complex  electric  field 
in  the  focal  plane,  specifically 


e>  cd 

F rj (js)  jj  /  <27> 

Once  the  approximated  Fourier  series  coefficients  are 
calculated  they  are  used  to  reconstruct  the  aperture  electric 
field. 

In  order  to  evaluate  the  validity  of  the  assumption 
that  the  modulus  of  the  focal  plane  field  can  be  approximate 
the  complex  focal  plane  field  in  equation  (26)  an  error 
analysis  is  shown.  The  expected  mean  square  error  (E(error)) 
between  the  focal  plane  electric  field  and  the  modulus  of  the 
focal  plane  field  is 
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£  [error]  =  [£  Jj/ 1  ~  Uf  (*(•!*)  j  (28) 

where 

U fU.Y')*  IliOb.ft'ljexrfjHn-M)  (29) 

and 


$(xf,yf)  =  phase  of  the  focal  plane  field 
The  expected  error  reduces  to 

E  [error]  =  [  f  jf  / U(  (*, ,/.)/  //' 

-2E [iif lUifa.p)!  | i-cos(M.*)}|<k«k  (30) 

The  phase  of  the  focal  plane  field  for  arbitrary  aperture 
fields  is  not  known.  In  order  to  gain  insight  into  this 
error  analysis  the  phase  of  the  focal  plane  field  is  modeled 
as  a  Random  Guassian  Process.  It  is  assumed  that  the  mean 
of  the  aperture  phase  is  zero  and  the  variance  is  sigma 
squared  ( ar  ) . 

E  [*(*<>/*)]  =  0  (3D 
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and 


E  j~  i^)~]  =  cr2,  (32) 

Note  that  the  expected  value  of  the  phase  term  (exp(- j^(xf ,y^) ) ) 
is  the  definition  of  the  characteristic  function,  and  its 
expected  value  is 

[  =exp(-i0l)  (33) 

The  expected  value  of  the  cosine  term  in  equation  (30)  is 

E[cos  (4(x,1y,))j  =  Re  [e xi°  (34) 

So,  the  expected  mean  square  error  between  the  focal  plane 
electric  field  and  its  modulus  from  this  model  (Random 
Gaussian  Process)  is 

09  UO 

fM =  <3  f/  I  uf  I 

I  I-  Re  [expi-i^Yl  I  c/Xfofyf  (35) 

There  is  no  reason  to  believe  that  the  variance  of  the  focal 
plane  phase  is  small  (actually  it  is  expected  to  be  large) 
the  expected  mean  square  error  for  this  model  is 
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£ [error]:  2  f  $  I  Uf  c/y, 


(36) 


-«  -a# 

From  this  error  analysis  the  Phase  Numeric  1  Technique  is 
not  expected  to  he  a  good  signal  processing  technique  for 
the  retrieval  of  the  aperture  phase  from  intensity  data  in 
the  focal  plane  (Refs  11  and  12). 

The  second  numerical  approach,  developed  from  the 
Phase  Analytic  Technique,  to  determine  the  spatial  (relative) 
phase  of  the  aperture  field  is  called  the  Phase  Numeric  2 
Technique.  This  technique  is  a  brute  force  method  using  an 
International  Mathematical  Statistical  Library  (IMSL)  sub¬ 
routine  ZXMIN  to  estimate  the  complex  Fourier  series 
coefficients  (Ref  13).  These  Fourier  series  coefficients 
are  then  used  to  reconstruct  the  aperture  field.  This 
technique  is  developed  in  this  paper  for  an  aperture  field 
which  is  a  function  of  one  spatial  variable  (one-dimensional 
fields).  This  analysis  can  be  expanded  for  two-dimensional 
fields . 

The  subroutine  ZXMIN  minimizes  a  user  supplied  function 
of  N  variables  using  a  quasi-Newton  method  (Ref  14) .  The 
user  supplied  subroutine  is  a  cost  function  which  first  uses 
the  N  variables  (c(N))  estimated  by  ZXMIN  to  form  Fourier 
series  coefficients  in  rectangular  (case  a)  and  polar  (case  b) 
coordinates,  specifically,  case  a 
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(37) 


pN  RECTx 


C(T)  ♦ 


and  case  b 


FM  P OLfiRx 


Ci i)  exp(jc^r+  £)) 


(38) 


These  Fourier  series  coefficients  are  then  used  to  approx¬ 
imate  the  focal  plane  wavefront 


(39) 


The  cost  function  is  then  generated  point  by  point 


J  -  £  j  (40) 

X 

rfetccr«< 

where 


1^(1)  =  measured  focal  plane  intensity  at  point  X 

If(I)  =  calculated  focal  plane  intensity  from  the 
approximated  Fourier  series  coefficients  at  point  I 


and 


£(*>• 


(41) 


The  IMSL  subroutine  ZXMIN  iterates  varying  the  N  variables 
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(Fourier  coefficients)  minimizing  the  cost  funct  (J),  as 
defined  in  equation  (40) .  For  a  more  thorough  explanation 
of  the  IMSL  subroutine  ZXMIN  refer  to  reference  14,  "Fortran 
Subroutines  for  Minimization  by  Quasi-Newton  Methods",  by 
R.  Fletcher. 

The  number  of  Fourier  series  coefficients  calculated 
are  N/2  since  these  coefficients  have  real  and  imaginary 
parts.  Once  the  optimized  Fourier  series  coefficients  are 
determined  the  aperture  wavefront  can  be  computed  by  inverse 
Fourier  transforming  the  optimized  focal  plane  field. 

This  technique  is  similar  to  one  tried  by  W.H.  Southwell 
(Ref  6),  where  Zernike  circle  polynomials  are  used  to  approx¬ 
imate  the  phase  of  the  aperture  field  instead  of  the  Fourier 
series  expansion  method  used  in  this  research  study. 

Autocorrelation  Function.  In  order  to  decrease  the 
number  of  variables  that  have  to  be  optimized  using  the 
Phase  Numeric  2  Technique,  the  autocorrelation  function 
of  the  aperture  field  was  investigated.  The  reasoning  is 
that  the  use  of  more  of  the  available  information  could 
result  in  the  reduction  of  computational  time  and  could 
improve  the  convergence  properties  of  the  Phase  Numeric  2 
Technique.  This  technique  is  developed  in  this  report  for 
aperture  fields  which  are  a  function  of  one  spatial  varia¬ 
tion  (one-dimensional  fields). 

The  autocorrelation  function  (R(Ax))  of  the  aperture 
field  is 
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(42) 


rm*z"{i  z  (u„m) rj 

The  Fourier  transform  of  the  aperture  field  squared 

'  p 

( (?(Ua(xa,ya) ) )  )  in  the  focal  plane  is  a  measured  quantity 
within  a  known  constant,  specifically 


R(*x)  =  7*'1  [  (x,))  ^ 

The  auto  correlation  theorem  states  that  if  g(x)  and  G(f) 
are  Fourier  transform  pairs  then 

7  f  J  5(rh*^'x)  ~  I G  <44) 

(Ref  9) .  Using  the  representation  for  the  aperture  wavefront 
in  equation  (12)  and  equation  (16)  then  substituting  in 
the  autocorrelation  theorem  (equation  (44))  we  find 


R(ax)  = 


(45) 


After  preforming  the  integration  in  equation  (45)  the  above 
autocorrelation  function  (equation  (45))  reduces  to 
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R  M «  Z.  /  Fh  iZe*r  (eL^r  Ax)  («) 

From  equation  (46)  the  magnitude  of  the  Fourier  series  coef¬ 
ficients  squared  be  found 

|F„r^wZ  Ri**)  ™ 

specifically 


MX- <“>] 


Equation  (48)  is  a  direct  method  to  compute  the  magnitude 
of  the  Fourier  series  coefficients  which  can  be  used  in 
the  Phase  Numeric  2  Technique.  Using  this  result  only 
N  phase  terms  have  to  be  optimized  for  N  Fourier  series 
coefficients . 


Gerchberg-Saxton  Algorithm 

The  Gerchberg-Saxton  Algorithm  is  an  iterative  method 
developed  to  solve  the  phase  retrieval  problem  in  electron 
microscopy  (Ref  7).  Adapted  to  solve  the  imaging  sensor 
phase  retrieval  problem  the  Gerchberg-Saxton  Algorithm 


(Fig  4)  starts  with  an  initial  guess  of  the  aperture  wave- 
front  (^ai(xa)).  The  approximated  focal  plane  wavefront 
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(F*i  (xf.) )  is  the  Fourier  transform  of  the  initial  guess  of 
the  aperture  wavefront 


F,  (Xf) 


(49) 


This  approximated  focal  plane  wavefront  is  then  modified  by- 
dividing  it  by  its  modulus  (normalizing  intensity)  and  then 
multiplying  this  function  by  the  measured  modulus  in  the 
focal  plane  ((Ira(xf))^) 


FMt)- 


f,  Qr<) 

I F,  W I 


Vz 


(50) 


This  new  focal  plane  wavefront  (F2(x^.))  has  the  modulus  of 
the  measured  focal  plane  wavefront  and  a  phase  which  is  a 
function  of  the  initial  guess  of  the  aperture  wavefront 
(fa^(xa)).  Now,  the  approximated  aperture  plane  field 
(fap (xg) )  is  the  inverse  Fourier  transform  of  the  modified 
focal  plane  field  (F2(xf)) 

t Uu-w/' (5D 

This  approximated  aperture  plane  field  (^a2^xa^ 
output  of  the  Gerchberg-Saxton  Algorithm  and  is  a  function 
of  the  initial  guess  of  the  aperture  wavefront  (^ai(xa)) 
and  the  intensity  measurement  (^(Xf))  in  "the  focal  plane 
of  the  imaging  sensor.  The  output  of  the  first  iteration 
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Figure  4.  Block  diai 


gram  of  the  Gerchberg-Saxton  Algorithm 


26 


(fa2(xa)),  is  now  used  as  the  new  guess  of  the  aperture 
wavefront  (fa^(x&)).  This  procedure  is  continued  until 
convergence  or  the  desired  accuracy  is  achieved. 

Error  Reduction  Technique.  J.R.  Fienup  modified  the 
Gerchberg-Saxton  Algorithm  with  an  input/output  technique 
which  is  called  the  "Error  Reduction  Technique”.  This 
technique  yields  a  solution  which  improves  with  each 
iteration,  reducing  the  mean  square  error  between  the 
focal  plane  intensities  (calculated  versus  measured),  thus 
the  name  Error  Reduction  Technique.  Instead  of  using  the 
output  estimate  of  the  aperture  field  (^a2^xa^  directly 
as  the  next  guess  of  the  aperture  plane  field  (£'  (x  )) 

3.1  3. 

and  continue  to  iterate  (Fig  5),  J.R.  Fienup  uses  fal(x&), 

A 

fa2 'xa  ’  otJect  constraints,  and  feedback  techniques  to 
compute  the  new  aperture  field  guess  and  then 

iterates  (Ref  8) . 

With  the  Error  Reduction  Technique  in  mind,  two 
feedback  techniques  are  postulated  that  use  object 
constraints  that  are  related  to  the  retrieval  of  the  phase 
of  the  aperture  field  from  intensity  measurements  in  the 
focal  plane  problem  that  is  being  addressed  by  this 
research  study  (Fig  5)*  The  object  constraints  that  are 
used  for  this  project  are  that  the  wavefront  being  approx¬ 
imated  must  have  a  modulus  of  one  across  the  aperture  and 
zero  elsewhere. 

The  first  feedback  technique  postulated  is  when  the 
output  aperture  field  (f&2 (x&) )  approximated  satisfies 


(fa2(xa))i  is  now  used  as  the  new  guess  of  the  aperture 
wavefront  ^al^xa^*  Th^s  procedure  is  continued  until 
convergence  or  the  desired  accuracy  is  achieved. 

Error  Reduction  Technique.  J.R.  Fienup  modified  the 
Gerchberg-Saxton  Algorithm  with  an  input/output  technique 
which  is  called  the  "Error  Reduction  Technique".  This 
technique  yields  a  solution  which  improves  with  each 
iteration,  reducing  the  mean  square  error  between  the 
focal  plane  intensities  (calculated  versus  measured),  thus 
the  name  Error  Reduction  Technique.  Instead  of  using  the 
output  estimate  of  the  aperture  field  (^a2^xa^  directly 
as  the  next  guess  of  the  aperture  plane  field  (x  )) 
and  continue  to  iterate  (Fig  5) »  J.R*  Fienup  uses  f  Jx), 

311  cl 

A  v 

fa2(xa)*  object  constraints,  and  feedback  techniques  to 
compute  the  new  aperture  field  guess  (fal(x&)),  and  'then 
iterates  (Ref  8). 

With  the  Error  Reduction  Technique  in  mind,  two 
feedback  techniques  are  postulated  that  use  object 
constraints  that  are  related  to  the  retrieval  of  the  phase 
of  the  aperture  field  from  intensity  measurements  in  the 
focal  plane  problem  that  is  being  addressed  by  this 
research  study  (Fig  5).  The  object  constraints  that  are 
used  for  this  project  are  that  the  wavefront  being  approx¬ 
imated  must  have  a  modulus  of  one  across  the  aperture  and 
zero  elsewhere. 

The  first  feedback  technique  postulated  is  when  the 
output  aperture  field  (T&2^xa^  approximated  satisfies 
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Figure  5.  Gerchberg-Saxton  Algorithm 
controlling  feedback 


the  object  constraints,  the  aperture  field  (fa^(xa))  pixel 
(realization  at  point  x  )  stays  the  same.  If  the  output 
aperture  field  approximation  (fgp (x&) )  does  not  satisfy 
the  object  constraints,  feedback  is  used  to  compute  the 
aperture  pixel  (fa^(xa))»  specifically 


~fo~l  (*«-') 

e/se 


/  f  fAj.(x't)  Satis  r/*S  Oftj e*.r 


(52) 


Beta  is  a  feedback  operator  and  is  defined  in  chapter 
four,  “Simulation  Results  and  Analysis". 

The  second  feedback  technique  postulated  differs  from 
the  first  feedback  technique  in  that  the  estimated  aperture 
wavefront  is  always  zero  outside  the  aperture.  When  the 
output  aperture  field  approximation  ^a2^xa)^  satisfies  the 
object  constraints,  the  aperture  pixel  (fai(xa))  stays  the 
same.  If  the  output  aperture  field  approximation 
does  not  satisfy  the  object  constraints,  feedback  is  used  to 
compute  the  aperture  pixel  (^ai^xa))*  The  addition  constraint 
is  that  if  the  aperture  field  outside  the  aperture  is  not 
zero  it  is  set  to  zero,  specifically 


O 

A 


■f<Ki  (*&.)  =  J  -Pa.1  (**) 

(**V  ;  else 


•  oorsioe  At>e&Ti>6e 

.if  fU  (jfa)  SATltfier 
j  oe>J £<-r  to*Jsr*Aintrs  (53) 
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The  best  feature  of  the  Gerchberg-Saxton  Algorithm  is 
that  the  Fourier  series  coefficients  are  not  computed, 
eliminating  the  errors  that  are  inevitable  whenever  a  series 
representation  is  artifically  truncated. 

In  this  chapter,  the  phase  retrieval  techniques 
investigated  in  this  research  project  were  developed.  The 
first  approach  undertaken  is  addressed  at  finding  an  analy¬ 
tical  solution  to  the  phase  retrieval  problem.  When  a 
closed  form  solution  could  not  be  found  to  this  approach 
two  numerical  techniques  were  developed.  The  Phase  Numeric 
1  Technique  used  analytical  results  from  the  Phase  Analytic 
Technique  and  the  assumption  that  the  focal  plane  field 
can  be  approximated  by  the  modulus  of  the  focal  plane  field 
in  order  to  retrieve  the  phase  of  the  aperture  wavefront. 

The  Phase  Numeric  2  Technique  also  used  analytical  results 
from  the  Phase  Analytic  Technique  and  a  parameter  optimiza¬ 
tion  subroutine  to  approximate  the  Fourier  series  coeffic¬ 
ients.  These  Fourier  series  coefficients  are  then  used  to 
infer  the  phase  of  the  aperture  field.  The  second  approach 
investigated  to  retrieve  the  phase  of  the  aperture  field  is 
based  on  the  Gerchberg-Saxton  Algorithm.  The  Gerchberg 
-Saxton  Algorithm  and  two  feedback  techniques  are 
developed  which  directly  retrieve  the  aperture  phase. 
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IV  Simulation  Results  and  Analysis 


The  previous  chapter  developed  various  techniques 
directed  at  solving  the  phase  retrieval  problem.  This 
chapter  defines  the  aperture  wavefronts  used  in  the 
simulations,  the  measure  of  performance,  and  presents 
the  data  and  performance  analysis  of  the  phase  retrieval 
techniques  that  are  simulated.  The  techniques  simulated 
are:  the  Phase  Numeric  1  Technique,  Phase  Numeric  2 
Technique,  and  the  Gerchberg-Saxton  Algorithm  with  and 
without  feedback. 

Aperture  Wavefronts 

The  phase  retrieval  techniques  studied  in  this 
project  were  simulated  as  a  function  of  one  spatial 
variation  (one-dimensional  fields) .  The  aperture  field 
consists  of  15  sampled  pixels  (across  the  aperture)  and 
the  focal  plane  consists  of  32  pixels.  The  aperture 
wavefronts  that  are  used  for  the  simulations  are  complex 
with  a  modulus  of  one 

CXp(jWNw)  (54) 

The  wavefront  functions  simulated  are:  a  constant  term 
(a  trivial  case),  a  linear  term  (tilt),  a  sine  term  (odd 
function),  a  cosine  term  (even  function)  and  a  random  term 
which  is  linearly  distributed  between  zero  and  one. 
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Figure  5  illustrates  the  relative  (spatial)  phase 
(w  (x  ) )  across  the  aperture  for  the  linear,  sine  and  cosine 
wavefronts.  The  trivial  constant  phase  term  is  not  shown. 
Figure  6  illustrates  the  randoml (x_ )  and  the  random2(x  ) 
distributions  used  in  the  simulations. 

Performance  Measures 

In  order  to  compare  the  performance  of  the  simulated 
phase  retrieval  techniques  the  mean  square  error  between 
the  phase  of  the  aperture  wavefront  and  the  estimated  wave- 
front  is  used 


Msr  ( /  uuo-  wfc(oi; 


where 


a.»C*Tl/*£ 


w(L)  =  phase  of  the  aperture  wavefront  at  pixel  L 

3. 

w  (L)  =  phase  of  the  approximated  aperture  wavefront 

cl 


at  pixel  L 


Now,  in  order  to  compare  the  relative  accuracy  of  the  sim¬ 
ulated  phase  retrieval  methods  with  various  aperture 
wavefronts,  a  normalized  mean  square  error  is  used 


N  M5f  = 


f  | 

Z  |w„CL)(l 


In  the  case  where  the  denominator  in  equation  (56)  is  equal 
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0*0 


to  zero,  the  denominator  is  set  to  one  (this  occurs  when  the 
phase  of  the  aperture  wavefront  is  exactly  zero  across  the 
aperture) . 

Phase  Numeric  Approaches 

The  phase  numeric  approaches  use  the  measured  intensity 
in  the  focal  plane  to  generate  the  Fourier  series  coeffi¬ 
cients  of  the  aperture  wavefront.  These  approximated  Fourier 
series  coefficients  are  then  used  to  retrieve  the  phase  of 
the  aperture  wavefront. 

Phase  Numeric  1  Technique.  The  Phase  Numeric  1  Tech¬ 
nique  makes  use  of  the  orthogonal  set  of  basis  functions 
found  in  the  Phase  Analytic  Approach  to  generate  the  Fourier 
series  coefficients.  The  Fourier  series  coefficients  that 
are  computed  from  simulating  this  technique  are  not  accept¬ 
able.  The  Phase  Numeric  1  Technique  worked  well  on  the 

A 

computation  of  the  FQ  coefficient,  but  the  higher  order 
Fourier  series  coefficients  computed  had  unsystematic  errors 
which  are  attributed  to  the  approximation  that  the  magnitude 
of  the  focal  plane  field  can  be  used  in  equation  (27)  for 
the  complex  focal  plane  field.  The  error  analysis,  on 
pages  17,  18,  and  19  .  indicate  that  the  above  approximation 
is  a  poor  one  and  the  simulations  prove  it  is  a  poor  one. 

Phase  Numeric  2  Technique.  The  Phase  Numeric  2 
Technique  is  a  brute  force  method  using  an  IMSL  subroutine 
(Ref  13)  to  optimize  the  complex  Fourier  series  coefficients 
by  minimizing  the  mean  square  error  between  the  focal  plane 
intensity  measured  and  the  focal  plane  intensity  calculated. 
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For  N/2  Fourier  series  coefficients  there  are  N  real  variables 
which  have  to  be  optimized,  since  the  Fourier  series  coeffi¬ 
cients  are  complex.  These  Fourier  series  coefficients  are 
calculated  in  rectangular  (FNRECT,  case  1),  and  then  in 
Polar  (FNPOLAR,  case  2)  coordinates.  The  number  of  Fourier 
series  coefficients  (#FN's)  are  varied  between  one  (Fq)  and 
nine  (FQ,  F±1»  F±2>  F+^f  and  F±^).  Table  1  lists  the  best 
(least  NMSE)  versus  the  worst  (largest  NMSE)  performance 
of  the  Phase  Numeric  2  Technique  for  the  Fourier  series 
coefficients  computed  simultaneously.  It  is  interesting 
to  note  that  the  worst  performance  was  achieved  for  the 
number  of  Fourier  series  coefficients  equal  to  three  or 
seven,  while  the  best  performance  was  achieved  with  the 
number  of  Fourier  series  coefficients  equal  to  five  (except 
for  one  in  the  trivial  constant  phase  case) .  There  does  not 
seem  to  be  much  difference  in  Normalized  Mean  Square  Error 
(NMSE)  or  the  iterations  (ITT)  required  when  the  Fourier 
series  coefficients  are  computed  in  polar  or  in  rectangular 
coordinates.  Note  that  computation  (iterations)  are 
stopped  when  the  approximated  Fourier  series  coefficients 
are  not  changed  (three  significant  figures)  by  the  Phase 
Numeric  2  Technique  subroutine  ZXMIN  for  two  consecutive 
iterations  or  on  the  1000  iteration,  whichever  occurs 
first.  The  average  execution  time  is  0.025cp  execution 
seconds  per  iteration,  and  with  the  routine  often  requiring 
the  full  1000  iterations  allowed  (and  still  not  totally 
converging  to  a  solution)  equates  to  250cp  execution 


Aperture 

FNRECT 

FNP0LAR 

Phase 

NMSE 

#FN's 

ITT 

NMSE 

#FN’s 

ITT 

1.  Const 

0.0 

1 

13 

0.0 

1 

13 

1.7E-2 

7 

1000 

6.0E-10 

7 

1000 

2 .  Linear 

2.2E-1 

5 

614 

2.3E-1 

5 

1000 

3.4 

7 

1000 

3.0E-1 

3 

465 

3.  Sine 

9.8E-1 

5 

1000 

9.8E-1 

5 

920 

1.6 

3 

250 

1.5 

3 

290 

4.  Cosine 

9.9E-1 

5 

530 

8.3E-1 

5 

1000 

1.0 

3 

224 

1.0 

7 

1000 

Table  1.  Representative  Simulation  Results  -  Best  vs.  Worst 
Performance  for  the  Phase  Numeric  2  Technique  (Simultaneously' 
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seconds. 


This  technique  did  a  fair  job  in  approximating  the 
linear  phase  aberation  (NMSE  =  0.2)  and  zero  phase  (NMSE  = 
IE-2)  but  had  unacceptable  performance  for  the  cosine  and 
the  sine  (NMSE  =1)  phase  aberations.  Increasing  the 
number  of  Fourier  series  coefficients  should  improve  the 
accuracy  achievable,  the  problem  is  that  as  the  number  of 
variables  that  have  to  be  optimized  is  increased  convergence 
is  more  difficult  to  achieve  (while  increasing  round  off  errors 
in  the  IMSL  subroutine  ZXMIN) . 

In  an  attempt  to  increase  the  number  of  Fourier  series 
coefficients  without  increasing  the  number  of  variables  being 
optimized,  at  any  one  time,  the  Fourier  series  coefficients 

A  A 

were  computed  serially,  first  FQ,  then  F+1  and  so  on.  Table 
2  lists  representative  simulation  results  showing  best  (least 
NMSE)  and  worst  (highest  NMSE)  for  the  aperture  wavefronts 
indicated.  It  is  interesting  to  note  that  the  NMSE  with 
nine  Fourier  series  coefficients  were  either  the  best  or  the 
worst  simulation  cases  for  seven  out  ot  the  eight  cases 
shown.  This  technique  of  computing  the  Fourier  series  coeffi¬ 
cients  serially  yielded  similar  results  to  the  technique 
where  the  Fourier  series  coefficients  were  computed  simul¬ 
taneously.  Note  again  that  optimization  is  stopped  when 
the  approximated  Fourier  series  coefficients  do  not  change 
(three  significant  figures)  for  two  consecutive  iterations 
or  on  the  1000  iteration  (per  set  of  coefficients  being 
optimized),  whichever  occurs  first.  This  technique  does 
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Aperture 

FNRECT 

FNPOLAR 

Phase 

NMSE 

#FN's 

ITT 

NMSE 

#FN's 

ITT 

1.  Const 

2.1E-37 

1 

13 

2.1E-27 

1 

13 

5.0E-10 

3 

447 

1.7E-18 

9 

643 

2 .  Linear 

4.3E-1 

9 

815 

3.8E-1 

9 

1064 

4.9E-1 

3 

321 

4.4E-1 

3 

258 

3.  Sine 

1.0 

1 

66 

1.0 

1 

59 

1.6 

9 

742 

1.6 

9 

984 

4.  Cosine 

1.0 

1 

75 

1.0 

1 

58 

1.0 

9 

768 

5.1 

9 

929 

Table  2.  Representative  Simulation 

Results  - 

Best  vs 

.  Worst 

Performance  for  the  Phase  Numeric  2  Technique  (Serially! 


* 


l  ? 
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a  fair  job  in  the  approximation  of  the  linear  phase  function 
and  a  good  job  of  approximating  the  trivial  zero  phase  test 
case.  This  technique  did  an  unacceptable  job  in  approximating 
the  sine  and  cosine  aberation  cases. 

The  problem  with  computing  the  coefficients  serially 
is  that  errors  are  compounded,  specifically  the  error  in 

A 

the  FQ  coefficient  causes  errors  in  the  optimization  of  the 

A 

F±1  coefficient  and  so  on.  The  data  indicates  that  if  the 
number  of  Fourier  series  coefficients  computed  is  increased 
beyond  the  nine  computed,  the  performance  of  the  linear 
aberation  phase  retrieval  will  be  increased,  since  the  best 
performance  was  achieved  when  the  maximum  number  of  Fourier 
coefficients  (nine)  were  used.  The.  average  execution  time 
is  0.025cp  execution  seconds  per  iteration. 

To  summarize,  the  Phase  Numeric  2  Technique  is  a  fair 
method  to  compute  linear  phase.  Figure  8  illustrates  and 
compares  the  best  cases  for  linear  phase  retrieval;  FNRECT 
(#FN's  are  five)  for  Fourier  series  coefficients  that  are 
computed  simultaneously  and  FNPOLAR  (#FN's  are  nine)  for 
Fourier  series  coefficients  that  are  computed  serially. 

The  NMSE  comparisons  and  Figure  8  indicate  that  it  is 
better  to  compute  the  Fourier  series  coefficients 
simultaneously,  rather  than  serially. 

Autocorrelation  Function.  Several  computer  simulation 
runs  were  made,  using  the  autocorrelation  function  to  approx¬ 
imate  the  magnitude  of  the  Fourier  series  coefficients  and 
the  Phase  Numeric  2  Technique  to  optimize  the  phase  of 
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these  Fourier  series  coefficients.  In  every  case  tried  the 

A 

magnitude  of  the  Fq  coefficient  is  exact,  while  the  higher 
order  Fourier  series  coefficients  had  errors  in  them. 

These  errors  caused  the  Phase  Numeric  2  Technique  to  have 

A 

convergence  problems .  The  NMSE  was  one  only  when  the  FQ 
coefficient  was  used  to  reconstruct  the  aperture  wavefront, 
and  always  greater  than  one  when  the  higher  order  Fourier 
series  coefficients  were  used  to  reconstruct  the  aperture 
wavefront. 

The  analytical  solution  using  the  autocorrelation 
function  is  exact.  The  problem  with  the  computation  of 
the  magnitude  of  the  Fourier  series  coefficients  is  in 
the  numerical  technique  used.  The  width  of  Ax  in  equation 
(47)  was  varied  from  1  to  1000  per  pixel  and  it  was  found 
that  the  higher  order  Fourier  series  coefficients  are 
very  sensitive  to  the  width  of  Ax,  indicating  that  this 
is  where  the  numerical  problem  is.  Additional  debugging 
was  stopped  because  of  lack  of  time. 

Gerchberg-Saxton  Algorithm 

The  data  from  the  Gerchberg-Saxton  Algorithm  simulations 
indicate  that  it  has  good  convergence  properties.  The 
Gerchberg-Saxton  Algorithm  (without  feedback)  had  the 
lowest  average  normalized  mean  square  error  (NMSE)  and 
the  fastest  convergence  of  the  techniques  simulated.  The 
Gerchberg-Saxton  Algorithm  with  the  first  feedback  technique 
had  poor  performance,  never  converging  to  a  solution  with 
the  normalized  mean  square  error  fluctuating  on  every 
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iteration.  The  Gerchberg-Saxton  Algorithm  with  the  second 
feedback  technique  had  a  low  average  normalized  mean  square 
error  and  did  not  converge  as  quickly  as  the  technique 
without  feedback. 

Gerchberg-Saxton  Algorithm  (without  feedback) .  The 
Gerchberg-Saxton  Algorithm  has  good  convergence  properties. 
Table  3  contains  the  simulation  results  of  the  Gerchberg 
-Saxton  Algorithm.  The  initial  guess  of  the  aperture  field 
(^a^(xa))  simulated  ares  (l)  1  +  jO,  and  (2)  1  +  jrandom(xa). 
The  distribution  of  the  random  function  is  shown  in  Figure  6. 
The  simulation  results  indicate  that  the  better  the  initial 
guess  of  the  aperture  field  is  the  better  the  normalized 
mean  square  error  (NMSE)  of  the  solution  is  and  the  faster 
the  algorithm  converges.  The  actual  performance  of  the 
algorithm  is  better  than  Table  3  indicates  in  certain 
cases  (note  l).  In  the  cases  of  the  Linear  and  the  Cosine 
aperture  phase  aberations  the  approximated  aperture  field 
has  bias  errors.  Figure  9  compares  the  phase  of  the  aperture 
field  with  the  Cosine  aperture  phase  aberation  to  the 
approximated  aperture  phase  using  the  Gerchberg-Saxton 
Algorithm,  the  bias  error  is  apparent.  These  bias  errors 
are  unimportant  since  this  research  is  interested  in  retrieving 
the  relative  (spatial)  phase  of  the  aperture  field.  In  these 
cases,  when  the  bias  errors  sire  removed,  the  normalized  mean 
square  error  is  at  least  one  order  of  magnitude  better  than 
indicated. 
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Aperture  Phase 

NMSE 

#ITT 

Init  guess  of 

Aperture  Field 

1. 

Const 

3.8E-29 

1 

1+jO 

9.5E-3 

1 

l+jrandoml (x_ )/10 

cl 

2. 

Linear 

1.7E-2 

1 

1+00 

1.8E-27 

3 

1+jO 

5.4E-2 

1 

1+jrandoml (x  )/l0 

3.7E-2 

5 

1+jrandoml (x~)/l0 

3.6E-2 

10 

l+jrandoml  (xT^/lO 

(note  l) 

3.4E-2 

50 

1+jrandoml  (x;~)/l0 

3. 

Sine 

7.0E-2 

l 

1+00 

1.6E-27 

3 

1+00 

1.0 

1 

1+0 randoml (x  )/l0 

9.1E-2 

10 

l+jrandoml  (xj~)/lO 

a 

4. 

Cosine 

1 

1 

1+00 

9.6E-2(* 

)  1 

1+  0 randoml ( x  )  /l  0 

(note  1) 

3.0E-2(* 

)  2 

1+O'randoml  (xT^/lO 

cL 

5. 

Randoml (x  )/l0 

8.8E-1 

1 

1+Orandom2 (x  )/lO 

1.6E-1 

8 

1+Orandom2 (x^j/lO 

d 

6. 

Randoml(x  ) 

8.0E-1 

1 

1+0  r andom2 ( x  ) /l 0 

a 

7.9E-1 

2 

1+0  rand  om2  (xfO/lO 

9.3E-1 

5 

l+;jrandom2  (x~3/l0 

9.8E-1 

10 

1+Orandom2  (xr3/l0 

d 

(*)  indicates  convergence 

to  the  conjugate 

(note  1)  indicates  bias 

errors , 

NMSE  is  at  least 

one  order  of  magnitude 

better 

Table  3.  Simulation  Results  -  Gerchberg-Saxton  Algorithm 
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In  the  cases  when  the  aperture  phase  aberation  is  an 
odd  function,  the  algorithm  converges  to  the  conjugated 
aperture  phase,  while  in  all  other  cases  the  algorithm 
converges  to  the  actual  phase  (not  conjugated) .  The 
literature  indicates  that  this  non-uniqueness  is  not  a 
problem  in  aperture  phase  retrieval  when  the  aperture  field 
is  a  function  of  two  spatial  variations  (Ref  8).  The 
simulations  performed  indicate  that  this  technique  converges 
very  quickly  for  the  first  few  iterations  and  then  slows 
down  for  the  remaining  iterations.  The  Linear  aperture 
phase  aberation  (Table  3»  case  2)  with  the  1  +  jrandoml(x  ) 
initial  guess  of  the  aperture  field  is  a  good  example  of 
this  characteristic;  after  the  first  iteration  (ITT)  the 
normalized  mean  square  error  is  5»4E-2,  after  the  fifth 
iteration  the  normalized  mean  square  error  is  3.7E-2,  and 
after  fifty  iterations  the  normalized  mean  square  error 
is  3.4E-2.  The  Gerchberg-Saxton  Algorithm  converges  to  a. 
good  solution,  with  a  normalized  mean  square  error  less  than 
1.6E-1,  in  less  than  10  iterations  for  all  cases  simulated 
except  for  the  Cosine  (case  4)  as  an  aperture  phase  abera¬ 
tion  with  the  initial  guess  of  the  aperture  field  of 
1  +  jO,  and  the  Randoml(x)  phase  aberation.  The  Gerchberg 
-Saxton  Algorithm  is  an  excellent  technique  to  estimate 
Linear,  Sine,  and  Cosine  aperture  phase  aberations,  and 
works  well  for  small  random  aperture  phase  functions 
(Table  3#  case  5)*  This  technique  requires  0.015cp 
execution  seconds  per  iteration.  If  it  is  estimated  that 
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that  10  iterations  are  required  for  convergence,  then  0.1 5cp 
execution  seconds  are  needed. 

Gerchberg-Saxton  Algorithm  with  Feedback.  The 
Gerchberg-Saxton  Algorithm  with  feedback  techniques  use 

object  constraints,  and  feedback  methods 
to  estimate  the  next  aperture  field  guess  (f'.(x  ))  as  a 
modification  to  the  Gerchberg-Saxton  Algorithm. 

The  first  feedback  technique  simulated  uses  feedback 
whenever  (inside  and  outside  toe  aperture)  the  object 
constraints  are  violated.  The  feedback  forms  simulated 
are 


P  (.Xck)  s 


;  F0f 

f  b(50O) 

;  foz 

]FZ  3 

f 

;  rei 

fb(j  <(*•.■>) 

The  output  aperture  field  is 

A 

■f«  M  *  s  (x*.)  <  6f«.)  (58) 

The  feedback  coefficient  (fb)  is  a  real  constant.  The 
initial  guess  of  the  aperture  field  is  either  1  +  jO  across 
the  aperture  and  zero  elsewhere,  or  1  +  jrandom(xQ)  across 
the  aperture  and  zero  elsewhere. 
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This  technique  would  iterate  to  a  fairly  low  normalized 
mean  square  error  and  then  on  the  next  iteration  fluctuate 
to  a  high  (above  one  in  some  cases)  normalized  mean  square 
error.  Table  4  lists  representative  convergence  charac¬ 
teristics  of  this  technique  for  the  cases  simulated.  In 
many  cases  simulated  the  lowest  mean  square  error  achieved 
was  on  the  first  iteration.  The  poor  convergence  charac¬ 
teristics  of  this  technique  make  it  unacceptable. 

The  second  feedback  technique  simulated  differs 
from  the  first  technique  in  that  the  guessed  aperture 
wavefront  is  always  zero  outside  the  aperture.  Simulations 
of  this  technique  indicate  that  it  has  good  convergence 
properties.  Table  5  lists  simulation  results  for  a 
feedback  function  where  the  imaginary  component  of  the 
output  aperture  field  (fa2(x&))  is  used  (equation  (57).  FB^) . 
This  feedback  technique  is  able  to  approximate  the  phase  of 
aperture  wavefronts  that  have  Linear,  Sine,  Cosine  and 
Constant  (trivial)  aperture  aberations,  with  a  normalized 
mean  square  error  that  is  less  than  6.1E-2  and  requiring 
about  10  iterations.  This  technique  could  not  retrieve 
the  random  phase  aberations  wells  for  the  case  simulated 
(Table  5,  case  5)  the  best  normalized  mean  square  error 
achieved  is  8.5E-1  on  the  fifth  iteration  (stopped 
converging  on  this  iteration) .  This  technique  requires 
0.015cp  execution  seconds  per  iteration. 

Interesting  characteristics  of  the  Gerchberg-Saxton 
Algorithm  with  the  second  feedback  technique  ares  that 
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#ITT 

Aperture  Phase 

Linear  (NMSE)  Sine  (NMSE) 

1 

5.4E-2 

9.6E-1 

2 

2.0E-1 

1.4E+0 

3 

1.7E-1 

3.4E-1 

4 

4.1E-2 

2.1E-1 

5 

6.5E-2 

2.4E-1 

6 

3.5E-2 

2.2E-1 

7 

7.3E-2 

1.5E-1 

8 

3.6E-2 

2.1E-1 

9 

5.8E-2 

1.9E-1 

10 

7.6E-2 

1 . 9E-1 

50 

7.3E-2 

1.8E-1 

100 

7.4E-2 

1 . 8E-1 

500 

8.0E-2 

1.6E-1 

Init  guess 

of  Aperture  Field  = 

1+ jrandoml (x  ) 

a 

Feedback  = 

B(^a2(xa))  =  1(FB2) 

Representative  Convergence  Characteristics  of  the 


Gerchberg-Saxton  Algorithm  with  the  First  Feedback  Technique 


Aperture  Phase 

NMSE 

#ITT 

Init  guess  of 
Aperture  Field 

1. 

Const 

3.8E-29 

1 

1+jO 

2.8E-1 

1 

1+ jrandoml (x  )/l0 

5.5E-3 

5 

1+ j randoml (xTO/lO 

2.4E-4 

20 

1+ jrandoml (xa)/lO 

1.4E-4 

30 

1+  jrandoml  (xa)/lO 

ci 

2. 

Linear 

1.7E-2 

l 

1+jO 

7.3E-3 

3 

1+jO 

1.2E-3 

10 

1+jO 

4 . 5E-4 

20 

1+jO 

3. 

Sine 

7.0E-2 

1 

1+jO 

1.0E-4 

2 

1+jO 

1.0E+0 

1 

1+ jrandoml (x  )/l0 

1.4E-1 

10 

1+ jrandoml (xa)/l0 

3. 

4. 

Cosine 

1.0E+0 

1 

1+jO 

9.6E-2 

1 

1+  j  rand  om 1 ( x  )  /l  0 

6.1E-2 

2 

1+j randoml (xa)/l0 

5. 

Randoml (x„)/l0 

8.8E-1 

1 

l+jrandoml(x  )/l0 

a 

8.5E-1 

5 

1+ jrandoml (xa)/l0 

8.5E-1 

10 

1+ jrandoml (xa)/l0 

SI 

.a 

Feedback  =  B(fa2 

(xa))  = 

1(FB3) 

Table 

*5.  Simulation  Results  - 

Gerchberc-Saxton  Alcorithm 

with  the  Second  Feedback  Technique 
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odd  aperture  aberations  (Table  5»  case  4)  converge  to  its 
unique  solution  and  retrieved  the  phase  of  the  aperture 
field  without  bias  errors,  unlike  the  Gerchberg-Saxton 
Algorithm  without  feedback. 

This  chapter  defines  the  aperture  wavefronts  used 
in  the  simulations,  the  measure  of  performance,  and  presents 
the  data  and  performance  analysis  of  the  phase  retrieval 
techniques  simulated.  The  Phase  Numeric  1  Technique 
yielded  unacceptable  performance.  The  Phase  Numeric 
2  Technique  is  able  to  approximate  linear  phase  aberations 
(tilt)  in  the  aperture  field,  but  could  not  approximate 
more  complicated  aberations.  This  technique  required 
approximately  250cp  execution  seconds  to  reconstruct  a 
linear  aberation.  The  Gerchberg-Saxton  Algorithm  yielded 
the  best  performance  of  the  techniques  investigated 
converging  to  a  solution  in  most  simulation  runs  within 
0.075cp  execution  seconds. 
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V  Conclusion 


Retrieval  of  the  phase  of  the  aperture  field  from 
intensity  measurements  in  the  focal  plane  of  an  imaging 
sensor  is  possible.  The  aperture  phase  although  not 
explicitly  present,  is  contained  in  the  intensity  dist¬ 
ribution  in  the  focal  plane  of  an  imaging  sensor.  Phase 
retrieval  techniques  studied  in  this  project  were  simulated 
for  an  aperture  field  which  is  a  function  of  one  spatial 
variation  (one -dimensional  fields).  Three  techniques 
simulated  are  able  to  retrieve  the  phase  of  an  aperture 
field  with  various  success;  these  techniques  are:  the 
Phase  Numeric  2  Technique,  the  Gerchberg-Saxton  Algorithm 
without  feedback,  and  the  Gerchberg-Saxton  Algorithm  with 
the  second  feedback  technique. 

The  Phase  Numeric  2  Technique  is  a  brute  force  method 
which  optimizes  the  Fourier  series  coefficients  as  a  means 
to  infer  the  aperture  phase.  This  technique  is  able  to 
retrieve  Linear  phase  aberations,  with  a  normalized  mean 
square  error  of  2.2E-1  (requiring  250cp  execution  seconds 
to  converge  to  a  solution) »  but  can  not  retrieve  more 
complicated  aperture  phase  aberations. 

The  Gerchberg-Saxton  Algorithm  (without  feedback)  has 
the  best  performance  of  the  techniques  simulated  in  this 
research  project.  The  Gerchberg-Saxton  Algorithm  converges 
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to  a  solution  with  a  normalized  mean  square  error  less  than 
1.6E-1  (in  less  than  0.15cp  execution  seconds)  for  all 
aperture  phase  alterations  simulated  except  for  aperture 
phase  alterations  with  large  random  fluctuations.  The 
Gerchberg-Saxton  Algorithm  with  the  second  feedback  method 
is  able  to  retrieve  the  aperture  phase  of  the  wavefronts 
simulated  with  a  normalized  mean  square  error  less  than 
6.1E-2  (in  less  than  0.15cp  execution  seconds),  except 
for  aperture  phase  aberations  with  random  fluctuations. 
Suggestions  for  Further  Study 

An  area  of  recommended  research  is  to  investigate 
the  technique  developed  in  this  study  for  the  computation 
of  the  magnitude  of  the  Fourier  series  coefficients  using 
the  autocorrelation  function  of  the  aperture  plane.  This 
information  from  the  focal  plane  intensity  measurements 
has  the  potential  of  improving  the  performance  of  the 
Phase  Numeric  2  Technique. 

Another  issue  worthy  of  investigation  is  to  determine 
how  the  Gerchberg-Saxton  Algorithm  Techniques  and  the 
Phase  Numeric  2  Technique  perform  with  an  aperture  field 
that  is  a  function  of  two  spatial  variations  (two-dimensional 
fields).  Additional  analysis  should  then  be  preformed  to 
determine  the  convergence  sentivities  of  these  approaches  in 
the  presence  of  noise  and  non-exact  intensity  measurements 
in  the  focal  plane  of  the  imaging  sensor. 
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